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Multifocal multiphoton microscopy (MMM) achieves fast imaging by simultaneously scanning multiple 
foci across different regions of specimen. The use of imaging detectors in MMM, such as CCD or CMOS, 
results in degradation of image signal-to-noise-ratio (SNR) due to the scattering of emitted photons. SNR 
can be partly recovered using multianode photomultiplier tubes (MAPMT). In this design, however, 
emission photons scattered to neighbor anodes are encoded by the foci scan location resulting in ghost 
images. The crosstalk between different anodes is currently measured a priori, which is cumbersome as it 
depends specimen properties. Here, we present the photon reassignment method for MMM, established 
based on the maximum likelihood (ML) estimation, for quantification of crosstalk between the anodes of 
MAPMT without a priori measurement. The method provides the reassignment of the photons generated by 
the ghost images to the original spatial location thus increases the SNR of the final reconstructed image. 

Multiphoton excitation fluorescence microscopy has inherent 3D resolution due to the nonlinear depend- 
ence of excitation upon the incident light distribution^'^. The excitation region is localized to a femtoliter 
volume at the focal point of a high numerical aperture objective lens. Multiphoton excitation fluor- 
escence microscopy is routinely used in a variety of tissue imaging applications due to its excellent imaging depth, 
high resolution, and lower photo-damage. However, one of the practical limitations of multiphoton excitation 
fluorescence microscopy is its imaging speed, typically up to a few frames per second. While this imaging speed is 
sufficient in many cases, several types of applications require faster systems. These applications include the 
measurement of dynamic processes including calcium signaling and action potential propagation^ high through- 
put image cytometric studies of tissue physiology'*'^, and clinical applications where the effects of physiological 
motion should be minimized. In vivo imaging of small structures, such as neuronal synapses, across a full 
neuronal arbor can also be problematic due to the lengthy anesthesia required to scan large volumes at high 
resolution^. 

High-speed multiphoton imaging has been previously implemented using several approaches. The first 
approach is based on using high-speed scanners such as polygonal mirrors^, resonant mirror scanners^ or 
acousto- optical deflectors (AODs)^"*^ instead of galvanometric mirror scanners used in conventional multi- 
photon microscopes. The high speed scanners can typically achieve a scanning speed of up to about 30 frames 
per second with a comparable imaging depth in tissues as conventional multiphoton microscopy. However, these 
high speed scanning systems require shorter pixel dwell time resulting in lower image contrast and poorer SNR. 
This can be partially compensated by increasing the excitation laser power, but it is not always possible due to 
photo-damage of the specimen and excitation saturation*^'*^. 

The second method is two-photon wide-field imaging*^'*^. With the temporal focusing method, a 3D resolved 
plane is generated at once by two-photon excitation instead of a focus eliminating the need for lateral scanner 
resulting in faster imaging speed. However, its performance is often limited by the lower axial resolution 
compared with conventional multiphoton microscopy and lower SNR due to the scattering of emission photons. 
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Figure 1 | Two experimental approaches for the estimation of a scattering matrix, (a) Sparsely distributed beads sample in tissue phantom (2% 
Intralipid emulsion), and scattering matrix measurement from its MMM image, (b) Single focus excitation on a real sample and scattering matrix 
measurement from its MMM image. 



The third method is multifocal multiphoton microscopy (MMM), 
which is an intermediate approach between the two previously 
described techniques With a lenslet array or diffractive optical 
element (DOE)^^'^'*, a number of foci are generated simultaneously 
and scanned together. The size of the whole scanning region corre- 
sponds to the sum of the sub-region scanned by each focus. 
Therefore, the imaging speed is improved proportionally to the num- 
ber of foci. One practical limitation for the imaging speed of MMM is 
available laser power. For typical titanium-sapphire laser with several 
Watts of output, approximately one hundred foci can be effectively 
generated for tissue imaging resulting in approximately two orders of 
magnitude improvement in imaging speed. 

Similar to two -photon wide-field imaging, the SNR of multifocal 
multiphoton microscopy is also limited by the tissue scattering of 
emission photons resulting in lower SNR compared with single focus 
scanning multiphoton microscopy. Simultaneous detection of mul- 
tiple foci requires detectors with spatial resolution that can distin- 
guish signals from all the foci at the same time. Most multifocal 
multiphoton microscopes use imaging detectors such as CCD or 
CMOS cameras. However, when emission photons generated at 
one focus are scattered in a turbid specimen and arrive at detector 
pixels that do not map to the focus location, the SNR of the image 
suffers. This lower SNR can be compensated by the use of multianode 
photomultipUer tubes (MAPMT) in a descanned detection geo- 
metry^^. The larger detection area of each anode greatly reduces 
the crosstalk between foci due to the scattering of emission photons 
and significantly improves the image SNR. 

However, with severe scattering in a turbid specimen, especially at 
larger imaging depth, scattered emission photons will still arrive at 
neighboring anodes, resulting in the formation of ghost images, 
which are duplicates of the image acquired by one focus visualized 
in neighbor sub -images. To remove the ghost images and to increase 
the SNR of the original image, these scattered emission photons 
should be reassigned to their original pixels. We have shown that 
this can be accomplished by estimating the scattering matrix at one 
location of the specimen that characterizes the crosstalk between the 
anodes of the MAPMT due to emission photon scattering. By apply- 
ing the inverse of the scattering matrix to the acquired image, an 
improved SNR picture free from ghost images can be produced. The 
estimation of the scattering matrix was done experimentally^^ (Fig. 1). 
The first method of the estimation is using a tissue phantom with 



similar optical parameter as real tissues containing sparsely distrib- 
uted fluorescence particles. These fluorescence particles are imaged 
by the MMM. The sparsity of the particles ensures that the ghost 
images of each particle can be identified. By measuring the relation- 
ship between particle intensities in the correct pixel and all its neigh- 
bors, the scattering matrix can be calculated. However, the phantom 
sample is often a poor model of a real tissue specimen on the micro- 
scopic level where there is significant heterogeneity, potentially lead- 
ing to a significant error in image post-processing. Second, the 
scattering matrix can be calculated with a real sample by exciting 
only one focus at a time in the same MMM system. This method 
gives the exact coefficients for the scattering matrix. However, since 
the scattering matrix is a function of specimen location, due to het- 
erogeneity, and depth, due to increased scattering, scattering matrix 
must be measured frequently. This labor intensive measurement pro- 
cess partly negates the speed advantage of applying MMM imaging. 

In this work, we demonstrate a novel image post-processing 
approach for MMM that allows estimation of the scattering matrix 
without any additional experimental measurement. The proposed 
approach is based on maximum likelihood estimation quantifying 
the crosstalk between the different anodes of the MAPMT based 
on the actual optical model of the MMM system. This approach 
reassigns the photons that are originally detected at wrong anodes, 
appearing as ghost images, to the correct anode location. Therefore, it 
greatly improves signal strength and simultaneously minimizes ghost 
images in the neighboring sub -images. Since the proposed approach 
uses only the images generated by the MMM without any tissue 
dependent a priori information, the characterization of the scattering 
matrix can be estimated at each location and depth. A mathematical 
model of the proposed approach is presented and validated first with 
simulation data where the results show the convergence of the 
method to the actual fluorophores structures. The performance of 
the proposed method is then employed for fluorescence beads and 
mouse neuron images where the results are quantitatively analyzed. 

Methods 

MMM configuration. We have developed two MMM systems; one with 45 |im foci 
separation, and the other with 85 [im foci separation. Fig 2 shows the schematic of the 
MMM with 85 [im foci separation based on a diffractive optical element (DOE) for 
multiple foci generation and the MAPMT in descanned detection geometry. The light 
source used is Chameleon Ultra II (Coherent, Santa Clara, CA). The excitation laser 
beam is expanded and illuminates the 8 X 8 or 4 X 4 DOE (customized, Holo/Or, 
Rehovoth, Israel) depending on the excitation wavelengths (8X8 for 780 nm and 4 X 
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Figure 2 | Schematic of MMM based on the MAPMT. 



4 for 910 nm). The beamlets are shrunk according to the size of x-y galvanometric 
mirror scanners (6215H, Cambridge Technology, Lexington, MA), and are 
overlapped on the scanning mirrors. The beamlets reflected by the scanners are 
expanded again to fill the back focal plane of a 20 X water immersion objective lens 
with 1.0 NA (W Plan-Apochromat, Zeiss, Thorn wood, NY), and enter the aperture 
with different entrance angles. The objective lens generates an array of 8 X 8 or 4 X 4 
excitation foci on the focal plane in a specimen. The image is formed by raster 
scanning of the excitation foci across the specimen using the scanning mirrors. The 
excitation foci are separated from each other by 85 |im, and each focus scans slightly 
larger than the area of 85 \im X 85 \im in the sample plane to facilitate the montage 
assembly. The scan area covered by the 8X8 foci is 680 |im X 680 \im with the 
85 |im separation, and 340 |im X 340 [im by the 4X4 foci. The emission photons are 
collected by the same objective lens and descanned by the scanning mirrors. The 
descanned emission beamlets are essentially stationary in relation to the scanning 
speed. The emission beamlets are reflected by a dichroic mirror (Chroma 
Technology, Bellows Falls, VT) toward the MAPMT and are focused onto the center 
of each corresponding anode of the MAPMT (H7546B-20, Hamamatsu, Bridgewater, 
NJ). An IR blocking filter (BG39, Chroma Technology, Bellows Falls, VT) and a short- 
pass filter (ET680sp-2p, Chroma Technology, Bellows Falls, VT) is installed before 
the MAPMT to block the excitation light. The MMM system with 45 \im foci 
separation has a very similar system configuration, but uses a different model Ti- 
Sapphire (Ti-Sa) laser (Tsunami, Spectra-Physics, Mountain View, CA) pumped by a 
continuous wave Nd:YV04 laser (Millennia, Spectra-Physics, Mountain View, CA) as 
a light source, a different 8X8 micro-lenslet array (1000-17-S-A, Adaptive Optics, 
Cambridge, MA), and a different 0.95 NA objective lens (XLUMPLFL20XW, 
Olympus, Melville, NY). In this case, the 8 X 8 foci with 45 |im foci separation cover a 
360 [im X 360 |im scan area in the specimen. 

Image reconstruction methodology. Let 0{x,y,z) be the fluorophores distribution in 
the specimen assuming constant irradiance per unit volume with the specimen. For 
two-photon point excitation at scan position x', the excitation intensity distribution at 
the specimen projected by an objective with a normalized 3D intensity PSF h(x) is 
E(x,x') = Eod(x — x')(S)h(x) = Eoh(x — x'). Here Eq is the maximum intensity and the 
symbol, (x), represents the 3D convolution operator. The fluorescence intensity 
generated at the specimen due to two-photon process is: 

F(x,x') = Elh\x-xyO(x) (1) 
The effect of emission photon scattering can be modeled generally by a emission PSF 
hm(x). Assume hm(x) = h(x) in the absence of scattering and h (x) represents the PSF 
for the scattered emission photon. 

hm(x) = (1 - 0L)h(x) + och'ix) (2) 

Without the emission scattering effect, the emission photon scattering strength, a, is 
zero. In generaP^ even in the presence of scattering, the amplitude of the modification 
term is generally small. Also, it has been observed that h'(x) typically has fuU-width- 
at-half maximum that is orders of magnitude larger than that of h(xy^. 



The intensity distribution at the image plane in epi-detection geometry is: 

J(x,)/,xO = [£o^'(^-x')-0(x)](x)/z,(f)|,=o (3) 

For the MMM system, let AT X AT equally spaced foci be arranged on a rectilinear grid 

with reference positions {x[^, ■^jv}> where ^ „ = {mA,nA,0), m and n are the 

foci indices, and A is the separation distance between foci. In order to generate an 
image using MMM, this grid of foci must be scanned to cover each A X A sub-regions 
for the corresponding axial plane. The signal for an MMM system at location {x,y) of 
an imaging sensor (CCD or CMOS) at any scan location {i,j,k) can be written as. 



lMMMix,y,i,j,k) = „ 



(4) 



: {[£o'/2^(^-^^,„-^y-0(% + 4)](x)/z,(x,^)}|. 



For a given z-plane, k, using an imager (such as a CCD or CMOS camera with sensor 
size M X M) that integrates the signal for all the lateral scan steps, the final single 
plane image can be written as: 



/ (x,y,k)=y^I {x,y,U,k) 

MMM ^ V ' / / J MMM ^ V ' V ' / 



(5) 



It is clear that in this case, the final image is "blurred" by the emission PSF and SNR is 
degraded. Instead of using an imager such as a CCD, we can also use a multianode 
PMT. By doing so there are two main differences^^: (i) the image is descanned in the 
multianode PMT to ensure that the foci are center at each anode of the PMT, and (ii) 
the signal at each anode is integrated at each scanning step. 



lMD{x,y,i,j,k)= ^ {[£o^^(^-^,„ -^//)'O(^ + 4)]®^m(^'^)}|z = 0,x = x + ;A/M,y=7+;A/M 
N,N J. 

= Y dx"dy"dz"Elh^{x" -mA-iA/m,y" -nA-jA/M,z") 



■ 0{x" ,y" ,Z" + kA')hrn {X - X" ,y - y" ,Z - z" ,k) | ,.o..=.+,a/m,,=,+;A/^ 
N,N . 

■■ J2 dx"dy"dz"Elh^{x" -mA-iA/M,y" -nA-jA/M,z") 

■ 0{x" ,y" ,z" + kA')hm(x + iA /M - x" ,y +)A /M -y\- z",k) 



(6) 
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For specimen with tissue-like scattering property, we may assume h'(x) to be broad 
and smooth on the length scale of A. Therefore, the integrand of the integral may be 
replaced by a set of constant values C(a,b),im,n),k that are related to the radially inte- 
grated PSF (Ho) of the scattered light. Here {a,b) are the anode index and {m,n) are the 
foci index. The elements of C matrix is depth dependent, specified by k for the given z- 
plane. Considering that the emission PSF can be further separated into a part taking 
into account of scattering, Eqn. (6) can now be rewritten as. 



(V) 



lMD(a,b,i,j,k) = (1 - a(k))ElHo | dx" dy" dz'h\x" ,y" ,z'yOix" +aA + iA/M,y" + bA +jA/M,z' + kA') 

N,N |. 

+ a(A:)£^ ^C(a,b),(m,n),k dx"dy"dz'h^{x" ,y" ,z'yO{x" + mA + iA/M,y" + nA+jA/M,z' + kA') 
Then equation (7) can be further simplified as: 



lMDia,b,i,j,k) = (1 - (x(k))HoI(a,b,i,j,k) + (x(k) ^ C(a,b),(m,n),k^(m,n,/,j,/c) (8) 



(i) Maximization of the likelihood function with respect to l{a,h,i,j,k) to get a 
better estimation for this variable: 



(ii) Maximization of likelihood function with respect to C to get a better estima- 
tion for this variable: 



Ul,Ho,C,oc,k) \ 



d^ls{I,Ho,C,a,k) 




C(a,b),(m,n),k = C, 



Note that we have generalized the emission photon scattering strength, a(/c), to be 
depth dependent but invariant over each 2D plane. The detection process consists of 
measured photon count by each anode {a,h) at each scan location (/, j, k), i.e. 
N{a,b,i,j,k). N{a,b,i,j,k) should follow Poisson statistics with a mean given by 
lMDici>ki,j,k). 

Maximum Likelihood Estimation. The log-likelihood function of the readout 
process, described in eqn. (8), can be written as. 

k= MaMj,k)lnlMDiaAiJ,k)-lMD(aAi,hk)] (9) 

a,b,i,j,k 

This log-likelihood function utilizes the actual optical model of the MMM for the 
maximum likelihood estimation of specimen fluorophore distribution. Our approach 
provides maximization of log-likelihood function and, thus, provides the 
quantification of the scattering matrix that is used for the reassignment of the image 
photons originally recorded as ghost image photons by the neighboring anodes of 
MAPMT to the correct spatial location. Therefore, photon reassignment can be 
casted as a blind deconvolution problem that seeks to recover the image lMDi^>hi,j,k), 
the depth dependent emission photon scattering strength, oc{k), and the mixing 
kernel, C(^a,b),{m,n),k given N{a,b,i,j,k). 

The original fluorophores distribution at the corresponding axial depth recorded 
by the MMM system can be considered as the 'first guess' for the iteration process and 
the maximization of log-likelihood function can be performed using a numerical 
optimization method which iteratively improves the strength of the actual signal by 
reassignment of ghost- image photons to the correct spatial location. We adopted 
Newton's method for maximizing the log-likelihood function and the iteration step 
for maximum likelihood estimation of fluorophores distribution can be written as: 



(iii) Maximization of likelihood function with respect to a to get a better estima- 
tion for this variable: 



] dh(I,Ho 



,C,a,k) \ 
\ L = 



j d^ls{I,Ho,C,a,k) \ 



(12) 



Here I^'^K C^^\ and a^^^ are the estimated fluorophores distribution, scattering 
matrix coefficients, emission photon scattering strength respectively at iteration 
step. The maximization process simultaneously updates the new estimate of 
C^^^^\ and a^^^^'' for iterative step k. The log-likelihood function is evaluated 
at each iteration step for maximum likelihood estimation of fluorophores distri- 
bution. The accuracy of this estimation partly depends on the initial estimates chosen 
and on the constrained parameters, i.e., positivity of the fluorophores. In the 
manuscript, we describe the process evaluating initial estimates of the unknown 
parameters from the recorded MMM image and show the applicability of the pro- 
posed method for simulation and experimental results. Fig. 3 provides the summary 
of the proposed photon reassignment algorithm for the MMM system. 

Results 

Simulation results. 4 Spots in a 2 by 2 MMM image. To test the 
feasibility of our approach, we started from the simplest case using 
simulations. We simulated a replica of a 4 bead image in a 2 by 2 
MMM [Fig. 4(a-(b))]. Fig. 4(a) shows the original image of 4 beads 



MMM image 
lMD(a,b, k) 



Cleaner Image 

(first guess) 
l^'\a,b,i,j,k) 


Scattering coefficients 
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Figure 3 | Summary of the photon reassignment process for MMM. 
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Figure 4 | (a) Original 4 beads image in a 2 by 2 MMM. (b) A scattering- affected image of (a); one sub-image contains one primary bead image and 
three ghost bead images with the specified proportions, (c) The processed image of (b). (d) Original alphabet image in a 2 by 2 MMM. (e) A 
scattering- affected image of (d) in the same fashion, (f) The processed image of (e). 



without ghost images, which would be the target of our simulation. 
For the MMM image, shown in Fig. 4(b), each sub-image contains 
one primary bead image with three ghost bead images, and for the 
sake of simplicity there is no overlap between the beads. Any 
scattering- simulated sub-image contained 100% of the intensity of 
the primary bead, 20% of the intensity of the adjoining beads, 
and 14% of the intensity of the diagonal bead. We added Poisson 
noise in the simulated image by using the Hmnoise function of 
MATLAB (MathWorks, Natick, MA). The proposed method 
requires the first guess of the cleaner image i^^\a,bj,j,k), scattering 

and emission photon scattering 



\a,h),{m,n),k^ 



matrix coefficients 

strengtha^^\ We chose the actual MMM image as the first guess 
for the cleaner image. In order to keep the proposed method 
adaptive, and to automatically generate the first guess for 
scattering matrix coefficients and emission photon scattering 
strength we used the following steps: 

(i) Identify the spatial position of the pixel (and its MAPMT 
anode) that contains the maximum intensity in the entire 
MMM image. 

(ii) Observe the intensities at the pixels corresponding to same 
spatial location of all the neighboring anodes identified in 
step (i). 

(iii) Generate the intensity ratios of the maximum value pixel to 
neighboring anode pixels. 

(iv) Repeat the steps (i)-(iii) for a few more values next to the 
maximum intensity, if possible, and take the mean of the all 
these values (coefficients of the scattering matrix). 

(v) More coefficients of the scattering matrix for the other anodes 
are calculated in the same manner described in the previous 
steps. These calculated coefficients are used for the first guess 
of the scattering matrix. In order to conserve the photon 
count, the sum of each row of scattering matrix (representing 
one anode under observation and its contribution in all 
neighboring anodes of MAPMT) is normalized to 1. 

(vi) The mean value observed in step (iv) is also considered as the 
first guess of emission photon scattering strength i.e.a^^^ 



We processed the simulated MMM image by initializing of the 
parameters as discussed above. The iteration process, as shown in 
Fig. 3, seeks to maximize the log-likelihood function and assigns the 
scattered emission photons of the ghost images to the original spatial 
locations. During the iterative process we defined two constraints; 
first total number of photons of the MMM images is constant before 
and after the processing, and second the positivity of the pixel values 
(as the fluorescence signal cannot be negative). The reconstruction 
process took around 10 iterations and Fig. 4(c) shows the processed 
image. It can be visually observed that all the ghost images are com- 
pletely suppressed in the final reconstructed image. The bead intens- 
ity is increased after processing, representing the reassignment of 
photons from the ghost image locations. Quantified results for the 
improvement in signal strength and signal-to-ghost image ratio 
(SGR) are presented later in this paper. 

Alphabets in a 2 by 2 MMM image. Next, we tested our reassignment 
method in a simulation of an alphabet image. The alphabet image 
contains certain shapes and overlaps that provide more complexity 
as compared to the simulated beads image. Again, an original image 
without ghost images was generated [Fig. 4(d)] as the reassignment 
target. Then, ghost images were imposed on all sub-images in the way 
that one sub-image in a 2 by 2 MMM image contains one primary 
alphabet image with the ghost images of the other three alphabets 
[Fig. 4(e)]. A given scattering- affected sub-image contains 100% of 
the intensity from a primary sub-image, 20% of the intensity from 
the two nearest neighbor sub-images along the horizontal and vertical 
directions, and 14% of the intensity of the farther neighbor along the 
diagonal direction. To start the iteration process the first guess of the 
scattering matrix and emission photon scattering strength was calcu- 
lated as before. The final processed image, shown in Fig. 4(f) shows a 
cleanly reconstructed alphabet similar to the original. Clearly, the 
effects of ghost images are significantly suppressed from the MMM 
images by using the proposed approach. As discussed before, the image 
intensity of the processed image is higher than the original simulated 
image because of the reassignment of the ghost image photons to their 
original spatial locations from which they were scattered. 
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Figure 5 | (a) Intensity comparison of the simulated beads MMM image (blue, main anode + ghosts, subjected to Poisson noise) and their processed 
images (red) showing photon conservation after the process, (b-d) Scattering matrix elements distribution at 2 by 2 MMM anodes, (b) Original 
distribution used in simulation, (c) estimated distribution of the processed bead image, and (d) the alphabet image. One sub-image contains primary 
element of scattering matrix from the corresponding anode of image, and three scattering matrix elements showing the leakage of emission photons to the 
neighbor anodes due to scattering. 



To evaluate the performance of the proposed method, we quanti- 
tatively analyzed the reconstruction results. Fig. 5(a) demonstrates 
photon conservation between the total intensity of the simulated 
MMM bead images and the processed images after photon-reassign- 
ment. The intensities of the simulated bead and its ghost images in 
the MMM image are summed and the sum is compared with the 
summed intensity of the same beads in the processed image. The 
error bars, in the MMM images, show the range of Poisson noise 
added in the simulation. Clearly, the reconstructed intensities at the 
different anode locations are within this Poisson noise range dem- 
onstrating that our algorithm conserves photons when signals from 
the ghost images are reassigned to correct spatial locations. Next, we 
compared elements of the scattering matrix corresponding to differ- 
ent anodes. The scattering matrix elements for the (1,1) anode should 
be [0.61 0.14 0.14 0.1] based on the preset parameters and the 
acquired emission photon scattering strength with normalization, 
and the same ratio follows for other anodes depending on their 
spatial locations. Our target is to compare the recovered scattering 
matrix elements based on the proposed approach with the original 
one used for simulation. Fig. 5(b) shows the original scattering 
matrix elements used for simulation of both bead and alphabet 
images. Fig. 5(c) and (d) show the reconstructed scattering matrix 
elements for the bead and alphabet images respectively. In both cases, 
the reconstructed scattering matrix matches very well the original 
values used to generate this simulation. 

The accuracy of the proposed method is evaluated with known 
underlying structures and with realistic noise level. We demonstrate 
that the proposed algorithm can accurately recover the expected 
features, and that the photon number is conserved between the 



simulated and processed images as expected for correct reassign- 
ment. The first guess of unknown parameters can often be estimated 
with good accuracy, as discussed above, resulting in a faster, more 
accurate MLE process. In order to further demonstrate the influence 
of the first guess we processed the MMM alphabet image using 
different first guess estimates and quantitatively compared the 
results. Fig. 6(a) shows the reconstruction results (same number of 
iterations used to process Fig. 4(f)) when using an incorrect first 
guess. For the incorrect first guess, we used half the intensity ratios 
used for the correct initial guess. Fig. 6(b) shows the quantitative 
comparison of the processed results from the incorrect first guess 
to those obtained using the correct first guess. It can be clearly 
observed that in the case of wrong first guess the signal recovery is 
not complete and photons corresponding to ghost images remain at 
the neighboring anodes. The recovery of the correct reconstruction 
results, for the wrong initial guess, can be possible by increasing the 
number of iteration. We observed that for the wrong guess (used for 
the Fig. 6(a)), the correct reconstruction result is achieved when the 
number of iterations is doubled. Therefore the effect of incorrect first 
guess mainly results in added computation time but not at the 
expense of an incorrect image. 

Experimental results. Fluorescent beads image in 6 by 6 MMM with 
85 fim foci separation. We tested our proposed approach with a 
fluorescent bead image taken from the 6 by 6 MMM (with foci 
separation 85 jim). The sample was imaged in the 8 by 8 MMM 
system described in Sec. 2.1. The image was trimmed to contain 
only 6 by 6 sub -fields for further image processing due to low 
signal levels at the edge of the field of view. The sample was 
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prepared with 10 and 15 |im diameter fluorescent latex microspheres 
(F8836 and F21010, Molecular Probes, Eugene, OR) immobilized in 
3D by a 2% agarose gel. A 2% fat emulsion (Microlipid, Nestle, Vevey, 
Switzerland) was added to mimic the scattering characteristic of a 
tissue specimen. 

The left column of Fig. 7 shows the acquired fluorescent beads 
MMM images at four different depths (0 jim, 40 jim, 85 jim, and 
135 |im). As expected, the scattering effect intensifies with increased 
imaging depth resulting in more prominent ghost images. To sup- 
press the effect of the ghost images and to reassign the photons of the 
ghost images to their original spatial locations, we used our proposed 
method. For the 6 by 6 MMM image, total 36 X 36 = 1296 coeffi- 
cients are necessary for the scattering matrix C. As the scattering 
effect is different with depth, it is necessary to separately estimate 
the first guess of the scattering matrix and emission photon scatter- 
ing strength for each depth. To start the photon reassignment itera- 
tion process, we used the same criteria described in section 3.1.1 to 
evaluate the first guess of the scattering matrix and scattering 
strength. The iteration process took around 30 iterations and the 
right column of Fig. 7 shows the processed images at their corres- 
ponding depths. Clearly, most of the strong ghost images were 
removed successfully, and only the primary bead images were left 
with better signal strength. To understand the performance of the 
proposed method on the signal levels of the imaging bead and its 
effect on the ghost image, we analyzed some of the relevant imaging 
parameters of the original and the processed MMM image. 

Fig. 8(a) shows the strength of the signal of the bead at its original 
spatial location before and after processing. As expected, the strength 
of the signal is higher for the processed images, particularly at the 
larger imaging depths since the processing provides the reassignment 
of photons from the ghost image locations to the original bead loca- 
tion and thus increases the strength of the signal. This improvement 
in signal is expected to be higher for the larger imaging depths 
because the higher scattering contributes to stronger ghost images. 



To analyze the performance of the proposed method more precisely, 
the ratio of the signal at original bead location to the signal detected at 
its strongest ghost image location was calculated, and termed signal- 
to-ghost ratio (SGR), as shown in Fig. 8(b). SGR was calculated for 
original and the processed images at each imaging depth. As 
expected, the SGR is highest for the shallowest imaging depth 
because of low scattering. The SGR of the original image, which is 
significantly lower for the large imaging depth because of severe 
scattering, is substantially recovered in the processed image. The 
increase of the SGR at 0 jam depth even with no intensity change 
is due to the suppression of the background noise during processing. 

Elements of scattering matrix at the bead image location anode 
and adjacent nearest neighbors were analyzed for different imaging 
depth, as shown in Fig. 8(c). As the imaging depth increases, due to 
increase in the emission scattering the scattering matrix elements 
become flatter so that the ratio of main anode to neighboring anodes 
decreases as previously investigated in Ref [25]. This effect is further 
elaborated by plotting the values of emission photon scattering 
strength acquired at different imaging depths, as shown in 
Fig. 8(d). As the imaging depth increases, the emission photon scat- 
tering strength increases exponentially as expected. 

Mouse brain imaging in 4 by 4 MMM with 85 /urn foci separation. We 
next tested our approach with a mouse brain image taken in the 4 by 4 
MMM system. Thyl-GFP transgenic mice^^ had surgery for cranial 
windows that were bilaterally implanted over the visual cortices 
between 6-8 weeks of age^^. These mice express green fluorescent 
protein (GFP) in a sparse pseudo-random subset of neocortical neu- 
rons. Imaging was performed on adult mice (>3 months) previously 
implanted with cranial windows. Mice were anesthetized with iso- 
flurane (3% for induction and 1.5% during imaging). Anesthesia was 
monitored by breathing rate and foot pinch reflex. The head was 
positioned in a custom made stereotaxic restraint affixed to a stage. 
Cell body and dendritic arbors of inhibitory neurons labeled with 



SCIENTIFIC REPORT: | 4 : 5 1 53 | DOI: 1 0. 1 038/srep05 1 53 



8 



(a) 



Original 



Processed 



UU 



200 



150 



100 

t 




0.2 



(b) 



Depth: 100(im 
^ Depth: 180|.im 
— Depth: IGO^m 




-1 0 1 

Distance from anode 




* a 

—Fitted 



50 100 150 200 
Depth (|im) 



Figure 9 | (a) A mouse brain image with the 4 by 4 MMM system. Left acquired images at 100 ^m and 190 ^m imaging depths, and Right are the 
corresponding processed images, (b) Plot of scattering matrix elements for an anode and its distribution at neighboring anodes corresponding to 
different imaging depths of mouse brain, (c) Plot of emission photon scattering strength as a function of imaging depth. 



GFP in layers 2/3 of visual cortex were imaged. The excitation wave- 
length was 910 nm, the laser power per focus was about 42 mW, the 
dwell time was 40 [is with 0.5 |im pixel, and the image size was 
340 jim X 340 jam with 4X4 foci. 

The left column of Fig. 9(a) shows the MMM images of the mouse 
brain at 100 |im and 190 jim imaging depths. Ghost images can be 
visually observed at both depths; they are, of course, more prominent 




Figure 10 | Plot of log-likelihood function for number of iterations. 



for the larger imaging depth. The proposed approach was used to 
suppress the effect of the ghost images and to reassign the photons of 
the ghost images to their original spatial locations. For the 4 by 4 
MMM image, total 16 X 16 = 256 coefficients are necessary for the 
scattering matrix C. We used the same criteria as described above to 
evaluate the first guess of the scattering matrix and emission photon 
scattering strength. Processed images are shown in Fig. 9(a), the right 
column at the corresponding depths. The accuracy of the approach 
can be seen in the removal of all ghost cell bodies. Further, the 
resultant final images contain realistic dendrite branches that are 
connected to the neuron body, and without unconnected branches, 
as would be expected from known neuronal anatomy. We chose this 
sparse image to test the recovery of real image features after reas- 
signment because we can easily trace the neuronal dendrite. With the 
expectation of similar behavior of imaging parameters as shown in 
Fig. 8, we analyzed the behavior of scattering matrix elements and 
scattering strength at different depths and Fig. 9(b) and (c) show the 
plot of their values as a function of depth. Similar to the result of the 
bead sample in section 3.2.1, the scattering matrix becomes flatter as 
the imaging depth increases showing more widely spread emission 
photons due to scattering. As expected, the emission photon scatter- 
ing strength fits well to an exponential function. 

Fig. 10 shows the plot of log-likelihood function versus number of 
iterations for the processing of mouse brain images as shown in 
Fig. 9(a). Convergence of the log-likelihood function starts at around 
15 iterations, which provides the final processed image. It should be 
noted that imaging time and image processing time are not equival- 
ent, especially for in vivo experiments. For example, the increase in 
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Figure 11 | Liver image with the 6 by 6 MMM system, (a) Left acquired images at 20 [im and 30 |^m imaging depths, and right are the corresponding 
processed images, (b) Intensity Une plots for original and processed images for imaging depths 20 [im and 30 [im respectively. 



imaging speed using MMM reduces the typical time, from approx. 90 
minutes (pixel residence time: 40 |is, 756 X 756 pixels scanned per 
XY plane using single focus excitation, 125 XY planes in Z stack) to 
approx. 6 minutes (16 pixels at the same residence time of 40 |is), for 
high resolution imaging of an entire neuronal volume in the brain of 
a living mouse. This reduction in image acquisition time is critical in 
terms of reducing anesthetic duration to minimize physiological 
influence and animal attrition. Further, we observed that with equal 
pixel residence time, the 4 by 4 MMM imaging system improves 
imaging speed by 16 fold compared to a typical single focus multi- 



Table 1 Contrast connparison of origina 


and processed liver 


images for imaging depths 20 |im and 30 


jim 


Original 


Processed 


Depth: 20 urn 0.62 


0.87 


Depth: 30 0.54 


0.81 



photon excitation fluorescence microscope, but maintains almost 
equivalent image SNR. The total computational time taken for pro- 
cessing two 4 by 4 MMM images was around 40 seconds using a 
standard desktop computer with i5 processor. This is about 10 times 
longer than the time for image acquisition but still substantially 
shorter than single focus acquisition. More importantly, the proces- 
sing time could easily be reduced to below image acquisition speed by 
parallelizing computation with graphical processing units (GPUs). 
With ever increasing and cheap computation power, we have no 
doubt that the bottleneck lies in physical image acquisition rather 
than image processing. 

Liver tissue imaging using 6 by 6 MMM with 85 fimfoci separation. 
We further tested our approach with liver tissue imaging containing 
significantly denser features. The imaging was performed with a 6 by 
6 MMM system with 85 |im foci separation. The livers of male 
Wistar rats with an average weight of 200 g that underwent bile duct 
ligation surgeries were harvested in accordance to the Institutional 
Animal Care and Use Committee regulations at the Biological 
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Resource Center of A* STAR Singapore^^. The liver tissue specimens 
were obtained from the paraffin preservation of the Hvers after fixa- 
tion by cardiac perfusion with 4% paraformaldehyde. 

Fig. 11(a) shows the original MMM image and the processed 
image of the liver tissue sample at 20 |im and 30 jam imaging depths. 
The ghost images in MMM are not so prominent. Because of the 
dense features scattering mainly generates a relatively uniform back- 
ground haze at the neighboring anodes. In this case, the proposed 
method suppresses the photons in the uniform background haze and 
reassigns them to the original anode location. Fig. 11(b) shows the 
intensity line plots for both imaging depths. It can be clearly observed 
that for the processed images the background noise is suppressed and 
image features are extracted with better signal strength. Contrast of 
the features, highlighted in Fig. 11(b) at the 20 jim and 30 jim 
imaging depths, is calculated for the original and processed image. 
The quantitative comparison of contrast is shown in Table 1. Clearly, 
the proposed method improves feature contrast; essentially because 
of the improvement in signal strength due to concurrent suppression 
of the background haze at the respective anode and more effective 
utilization of the available image photons. 

Fluorescent beads and Mouse brain imaging in MMM with 45 fimfoci 
separation. We tested our proposed approach with a fluorescent 
beads taken from MMM system with 45 jam foci separation. The 
reduced separation of the foci, for the same scattering conditions 
increases the crosstalk between foci and results in stronger ghost 
images. The bead sample was prepared with 4 |im diameter fluor- 
escent latex microspheres (F8858, Molecular Probes, Eugene, OR) 
immobilized in 3D by 2% agarose gel. A 2% intralipid emulsion 
(Liposyn III, Abbott Laboratories, North Chicago, IL) was added 
to mimic the scattering characteristic of a tissue specimen. 

Fig. 12(a) shows the raw fluorescence bead image taken from the 8 
by 8 MMM with foci separation of 45 |im at images at a depth of 
150 |im. This image was trimmed to contain 6 by 6 sub-fields for 
further image processing. Fig. 12(b) shows the processed image after 
photon reassignment. As expected, due to the smaller foci separation 
emission photons were scattered more into the adjacent neighboring 
anodes of the MAPMT resulting in more pronounced ghost images. 
We used our proposed approach to suppress the ghost images in such 
a severe scattering case. From the processed image shown in 
Fig. 12(b), we can visually observe that the ghost images are signifi- 
cantly suppressed, and the remaining dim signal represents beads 
situated in neighboring axial planes. Despite the same sample scat- 
tering due to the smaller foci separation, emission photon scattering 
strengths were 0.98 for 45 |im foci separation at 150 jim imaging 
depth as compared to 0.88 for 85 |im foci separation at 135 |im 
imaging depth. 

Mouse brain imaging in MMM with 45 /nmfoci separation. Finally, 
we tested our approach for a mouse brain image taken from the 
MMM system with 45 |im foci separation. For the mouse brain 
imaging, a Thyl-GFP transgenic mouse^^ was deeply anesthetized 



with 2.5% Avertin (0.025 ml/g i.p.) and transcardially perfused with 
PBS, followed by 4% paraformaldehyde. Its brain was dissected and 
placed overnight in cold 4% paraformaldehyde. 300 |im thick cor- 
onal sections were sectioned with a vibrotome, then mounted and 
coverslipped on microscope sUdes using adhesive siUcone isolators 
(JTR20-A2-1.0, Grace Bio-Labs, Bend, OR). 

Fig. 13(a) shows the mouse brain image taken at 90 |im imaging 
depth using the 8 by 8 MMM system with 45 |im foci separation. 
This image was trimmed to contain 6 by 6 sub -fields for further 
image processing. It can be observed that the ghost images from 
neighboring foci are substantial and these ghost images are intri- 
cately mixed with the "real" image features. Furthermore, a uniform 
background haze due the dense features of sample, as discussed 
before, can also be clearly visualized. We processed the image by 
using the same convergence criteria and the same method for estim- 
ating the first guess of the scattering matrix and scattering strength 
parameters. The processed image is shown in Fig. 13(b). Apparently, 
the ghost images and background haze are simultaneously sup- 
pressed from the processed image while the photons are assigned 
back to the real image. Thus it provides better contrast and SGR. 
Fig. 13(c) and (d) show the intensity line plots in the x- and y- 
directions, along the dotted lines shown on Fig. 13(a). 
Improvement of the signal and suppression of the ghost images 
can be clearly observed from the intensity line plots. One particular 
region highlighted in the images shows the overlap of ghost image 
with true image features. The recovery of the image features while 
simultaneously suppressing the ghost image can be observed in the 
processed image. We further show other imaging parameters e.g. 
strength of the neuronal signal at its original spatial location before 
and after processing (Fig. 13(e)), and distribution of the neuronal 
scattering matrix at the correct image anode location and as a func- 
tion of the distance from this location (Fig. 13(f)). 

Conclusions 

MMM has the advantage of improved imaging speed critical for in 
vivo imaging applications. However, due to scattering in turbid speci- 
mens MMM images suffer from ghost images especially at higher 
imaging depths. Previously, these ghost images have been processed 
using the scattering coefficients acquired experimentally. However, 
the experimental estimation requires frequent measurement with 
different specimens, different locations, and at different imaging 
depths. In this paper, we propose estimating the scattering coeffi- 
cients for the photon reassignment process based on the maximum 
likelihood (ML) estimation. This image post-processing method 
increases the SNR and SGR of the final processed image by reassign- 
ing the scattered photons to the original spatial locations, and also 
avoids experimental calibration of the scattering matrix resulting in 
greater experimental efficiency. It is also important to note that most 
experimental methods for calibrating the scattering matrix can only 
be performed on tissue phantoms or at a few locations in the spe- 
cimen, so that inevitable errors in matrix estimation in other loca- 
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Figure 13 | (a) Image of GFP expressing neurons in a mouse brain acquired using a 6 by 6 MMM system at 90 |am depth, and (b) is the corresponding 
processed image, (c) and (d) Intensity line plots in x- and y- directions for the mouse brain image, (e) Signal and SGR comparison for original and 
corresponding processed image, and (f) plot of the scattering matrix elements for an anode and its distribution at neighboring anodes. 



tions will greatly degrade reassignment accuracy and introduce sys- 
tematic error. This problem is completely avoided using our new 
approach where the only a priori parameters used are the physical 
properties of the optical system such as the objective numerical aper- 
ture and the light wavelengths. We have shown the feasibility of our 
approach with simulated results and confirmed that the algorithm 
can estimate the scattering coefficients and reassign ghost images to 
their correct location. We validated our approach with fluorescent 
bead samples in the turbid medium, as well as with in vivo mouse 
brain images. The processed images demonstrate that the scattered 
emission photons are reassigned to their original spatial locations 
resulting in higher signal. The SGR is improved by up to a factor of 
10. The algorithm also confirms that shorter foci separation gener- 
ates more crosstalk, and the emission photon scattering strength 
increases as a function of imaging depth. To minimize the crosstalk 



between foci and the resulting ghost images inherently, it is recom- 
mended to adopt a large field of view objective lens and maximize 
foci separation. However, in the case that only small field of view 
objective lenses are available or the region of interest is small, our 
proposed approach can be a good solution for image post processing 
and image recovery. It is also interesting that our approach can 
recover the scattering matrix and the emission photon scattering 
strength at each image location allowing us to quantitative recover 
tissue optical properties locally. 
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